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In this paper we shall propose a simple scheme for calculating Green's functions for photons 
propagating in complex structured dielectrics or other photonic systems. The method is based on 
an extension of the finite differMice time domain (FDTD) method, originally proposed by Yeed, also 
known as the Order-N methodu, which has recently become a popular way of calculating photonic 
band structures. We give a new, transparent derivation of the Order-N method which, in turn, 
enables us to give a simple yet rigorous derivation of the criterion for numerical stability as well as 
statements of charge and energy conservation which are exact even on the discrete lattice. We imple- 
ment this using a general, non-orthogonal co-ordinate system without incurring the computational 
ON . overheads normally associated with non-orthogonal FDTD. 

We present results for local densities of states calculated using this method for a number of sys- 
tems. Firstly, we consider a simple one dimensional dielectric multilayer, identifying the suppression 
in the state density caused by the photonic band gap and then observing the effect of introducing 
a defect layer into the periodic structure. Secondly, we tackle a more realistic example by treating 
a defect in a crystal of dielectric spheres on a diamond lattice. This could have application to the 
design of super-efficient laser devices utilising defects in photonic crystals as laser cavities. 
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One of the principle driving forces behind the recent flurry of research into photonic band gap materialstffl has 
been the potential for manipulating the spontaneous emission of an atom placed in a cavity in such a material. As 
was pointed out some years agoO a material with a periodically structured dielectric function (a photonic crystal) can 
have a profound effect on the density of states for photons within the material, in some cases leading to frequency 
I i windows for which no allowed photon states exist. These windows, or photonic band gaps, radically alter the emission 
properties of atoms. An excited atom that wants to emit a photon of a frequency within the band gap cannot do 
so — the photon has no states into which it can go and so is forced to from a new kind of atom-photon bound state. 
Even zero-point fluctuations are forbidden within the band gap. This has immediate implications for device physics. 
Placing an active device, such as a semiconductor laser, within a cavity in a photonic crystal offers the possibility to 
control unwanted spontaneous emission and allow emission only into the lasing mode, thus dramatically improving 
the efficiency of the devicaJ. 

Green's functions have the potential to play a central role in the theoretical investigations of these photonic systems. 
Not only are they a natural way to express key quantities such as the density of states, they are also are easily calculated 
within the framework of time domain methods such as the Order N technique that we shall be exploring in this paper. 
In systems where dissipation is present, such as those containing metals or lossy dielectrics, the Green's function 
approach is the only route to calculate quantities of physical interest. 

For the theorist, the challenge is to solve Maxwell's equations for such systems. Much progress has been made in 
the past few wars and several well established techniques have emerged. Probably the most widely used is the Plane 
Wave methocon. Simply put, this method involves expanding the electromagnetic fields as a sum of plane waves and 
recasting Maxwell's equations into the form of a eigenvalue problem to find the allowed eigenfrequencies. Though 
simple to implement and use, this method has the draw back that the time taken for the calculation scales as the 
cube of the number of plane waves used, so for complicated problems which require many plane waves this is a severe 
limitation. Another limitation is to systems whose dielectric functions do not disperse with frequency. Hence metallic 
systems are beyond the scope of this method. 

A second popular method involves working at a fixed frequency but instead of expanding the wavefield on a lattice 
in reciprocal space the wavefield is represented on the points of a real space latticeEJO. The resulting equations can 
be rearranged into the form of a Transfer Matrix which relates the fields in one layer of the lattice to the fields in the 
next. This method has proved extremely useful especially for systems involving metals where the dielectric constant 
is a function of frequency. Also, because the form the transfer matrix takes, connecting the fields on one surface to 
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the fields on another, calculations based on this method scale as the square of the number of real space points, rather 
than the cube. 

This may be an improvement over the.-plane wave method but is still worse than the optimal linear scaling with 
system size. However, it has been showno that, because Maxwell's equations are local then by working in the time 
domain instead of the frequency domain it is possible to obtain methods which scale as 'order iV' where N is the 
system size. We shall work with such a method in this paper in order to calculate photonic Green's functions and 
from those other quantities of physical interest such as densities of states. Previous work has_successfully calculated 
spontaneous emission rates from the density of states using a real space lattice formulationEj. However, that work 
was based in the frequency domain. By working in the time domain we are able to exploit the more favourable scaling 
law and consider larger, more complex systems such as photonic crystals with defects. 

II. METHODS 

The theory behind the methods used in this paper is somewhat similar in principle to the Finite Difference Time 
Domain method first introduced to the electrical engineering community by Yeeo in 1966, and first applied to the 
problem of photonic band structures by Chan et aln in 1995. Here we present a systematic derivation of the finite 
difference equations within a transparent formalism. The advantages of this new formalism over the traditional Yee 
approach are many. First, it makes clear how to find quantities which are exactly conserved even by the discrete 
equations. Hence we can identify the analogies to charge, energy density and the Poynting vector etc. Second, our 
formalism enables the precise analysis of the stability of the discrete equations to be made from a simple consideration 
of the approximations involved, sidestepping the usual, rather involved, Courant stability analysis. Third, and perhaps 
most importantly, the new formalism shows, for the first time, how to present the finite difference equations in a 
completely generaLco-ordinate system without incurring the computational overheads normally associated with non- 
orthogonal FDTDGJ. Finally, it is hoped that this transparency will make it simpler to extent the method to areas 
to which it has not yet been applied. We begin from the usual Maxwell's equations, neglecting any free charges or 
currents; 

<9E OH. 

VxH = e £(r)— V x E = - Mo /i(r)— (1) 

which on Fourier transforming into (K,w) space can be written, 

ik x H = — ieoeuj'E ikxE = +z/io^o;H (2) 

Next we wish to place this equations onto a discrete, real space lattice by replacing the derivatives with finite 
differences. We have some freedom here but must be careful — not all differencing schemes lead to stable equations. 
We wilL-take our lead therefore from the differencing scheme which has proved so successful in the transfer matrix 
methodO. We do this by introducing the following approximations to ui and k. For the terms which involve the 
electric field we use, 

(ifc x a) _ -i J-iujSt) _ i 

(3) 



la —i8t 
with similar expressions for and For the magnetic field terms we use, 

I _ g(-«fc x a) ^ _ e (iuiSt) 
k x ^>k~ = ; LO I ► LO — (4) 

la —idt 

etc. So making these approximations in equation (|^) we get, 

ikT x H = — ze £o; + E ik + x E = +i(i {nu~H (5) 

On Fourier transforming back into the (r, t) domain it is clear that these approximations are equivalent to taking a 
forward finite difference, A + , in place of derivatives of the electric field, and a backwards difference, A - , for derivatives 
of the magnetic field. 

A+F(r) = [.F(r + a)-.F(r)]/a (6) 
A-^(r) = [^(r)-^(r-a)]/a (7) 
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Putting it all together the discrete form of Maxwell's equations become, 

V+ x E(r, t) = -// M(r)A t -H(r, t) (8) 

V - x H(r, t) = e Q e(r) A+E(r, t) (9) 

where, 

-A+ A+ 

V+x = | A+ Z -A+ | (10) 
-A+ A+ 



(11) 

The approximations outlined in the previous paragraph place Maxwell's equations onto a discrete lattice of points 
which is uniform and Cartesian. However, this is not always convenient for the problems we may want to consider. It 
may be to our advantage to work in a co-ordinate system which is non-uniform or even non-orthogonal. Fortunately, 
as shown in detail elsewheraiil there is a simple result which allows us to map a completely arbitrary co-ordinate 
system onto a uniform Cartesian one as long as we introduce new, renormalised versions of the permittivity and 
permeability. In the generalised co-ordinates Maxwell's equations become, 

Z^M. eoi(r) ^5M (12) 





( o 













-A; 






A" 






n = -MoMr) 77 (13) 

Qo ot 

where, 

A+T{v,t)=T(r + Q 1 vi 1 ,t)-T{v,t) A+JT( r , t) = F(r, t + St) - F{r, t) (14) 



A~J : (r,t)=^(r,t)-J r (r-Q 1 u 1 ,t) A^( r , t) = F(r, t) - F{r, t - St) (15) 



etc, 



and 



&(t) = e(r)g^\m ■ u 2 x u 3 |^^i A « (r) = M ( r )^| Ul • u 2 x u 3 |^gi (16) 

*«(r) = a(r)^| Ul • u 2 x u 3 |^^i *« (r) = ^(r^m ■ u 2 x u 3 |^§i (17) 

The generalised co-ordinate system is defined by the three unit vectors, Ui , u 2 and u 3 which point along the generalised 
co-ordinate axis and the three lattice spacings Qi, Qi and Q% which define the spacing between discrete lattice points 
in each direction. In general, the Uj's and Q^s will themselves by functions of position. The metric tensor g lJ is 
defined such that (g^ 1 ) 1 -* = Ui ■ uy. As well as the permittivity and permeability being renormalised by the co-ordinate 
transformation, the fields themselves are also rescaled by the lattice spacings, 

Ei = QiEi Hi = QiHi (18) 

In order to obtain the equations that will link the fields at one time step to the fields at the next we introduce, 

H' = -^-H (19) 
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and let, 

a = QiUi b = Q 2 u 2 c = Q3U3 

Then after some rearranging we arrive at, 

Ei(r,t + St) = £i(r,t) + [e- 1 ^)] 11 [#£(r,t) - i%(r - b,t) - ff£(r,i) + tf 2 (r - c,t) 

+ [r 1 (r)] 12 (r, t) - H[ (r - c, t) - H' 3 (r, i) + i% (r - a, t)' 

+ [e-^r)] 13 f^(r,t)-^(r-a,t)-^(r,t)+^(r-b,t) 



£ 2 (r, t + Jt) = £ 2 (r, i) + [e-^r)] 21 [tf 3 (r, t) - H' 3 (r - b, t) - tf 2 (r, t) + H' 2 {v - c, t) 
+ [i~ 1 (r)] 22 (r, t) - H[ (r - c, t) - H' 3 (r, i) + H' 3 (r - a, t)' 
+ [^(r)] 23 [tf 2 (r, t) - H' 2 {v - a, i) - H[ (r, t) + H[(r - b, t) 

£b(r, t + St) =E 3 (r,t)+[e- 1 (r)] 31 [h' 3 (r, t) - H' 3 (r - b, t) - H' 2 (r, t) + H' 2 (r - c, t) 
+ [i~ 1 (r)] 32 (r, t) - H[ (r - c, t) - H' 3 (r , t) + H' 3 (r - a, t) 
+ [i~ 1 (r)] 33 i? 2 (r, t) - H' 2 (r - a, i) - ffj (r, t) + #J (r - b, t) 



(20) 



(21) 



(22) 



(23) 



#i(r,t + 5i)=irj(r,t) 



5< Cq 



Qo 



[fi- 1 (r)] 11 E 3 (r + b, t) - E 3 (r , t) - E 2 (r + c, t) + E 2 (r , t) 



WrM [A^'W] 12 Bi(r + c,t)- J B 1 (r,t)-£;3(r + a,t) + £!3(r,t) 



(5i Co 

~~Qo~ 



[A" 1 (r)] 13 £ 2 (r + a, t) - E 2 (r, t) - Eh. (r + b, t) + Eh (r, t) 



(24) 



J%(r,t + <fi)=#£(r,t) 



^r) [A" 1 W] 21 (r + b, t) - E 3 (r , t) - E 2 (r + c, t) + E 2 (r,t) 
Wo / L 

TT^ [A _1 ( r )] 22 [^i(r + c, t) - ^(r, t) - 4(r + a, t) + E 3 (v, t) 
wo / 1 

^r) [A" 1 (r)] 23 [4 (r + a, t) - E 2 (r, t) - E, (r + b, t) + E, (r, t) 
Qo ) L 



(25) 



H' 3 (r,t + St) = H' 3 (r,t) 



^r) [A _1 (r)] 31 [£ 3 (r + b, t) - 4(r, t) - E 2 (r + c, i) + E 2 (r, t) 
Wo J L 

TrO [A _1 ( r )] ^ N r + c '*) - (r, t) - ^ 3 (r + a, i) + £ 3 (r, i)" 

^) 2 [A-'W] 33 U 2 (r + a, t) - E 2 (r, t) - E^r + b, t) + E x (r, t) 
Wo J L 



(26) 
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These equations allow us to take an arbitrary set of electric and magnetic fields at some initial time t = and, subject 
to appropriate boundary conditions, calculate the fields at all later times. 

Notice that by incorporating all of the details of the generalised co-ordinate system into the new definitions of 
the permittivity and permeability we have, in effect, eliminated the extra computational overhead of the general co- 
ordinates. If the co-ordinate system changed from one place to another, and there is no reason why it should not, we 
would have to store the metric tensor at each point on the lattice. But since we have to store e and [i tensors at each 
point anyway, by combining the metric into our new i and ji we cut the amount of storage needed. Similarly, we also 
reduce the number of calculations required at each time step as we no longer need to worry about converting between 
covariant and contravariant vectors at every time step - this is all taken care of by e and fi. Hence, in contrast with 
previous methods, our non-orthogonal FDTD has no additional computational overhead compared to an orthogonal 
one, except for the initial setting up of e and fx. 



A. Stability Criterion 



These equations give a stable updating procedure for the fields if the time step is kept sufficiently small. The 
criterion is easy to find. Starting from the approximations we made for k and u>, (equations || and |]) the free space 
dispersion relation, to 1 



■ 2 — Cgk 2 , becomes, 



St 2 



sin 



Lo5t 



4c 2 , 



1 



Qx k a 



1 



1 



Q3 k z 



The condition that the maximum value of the right hand side must correspond to a real frequency gives, 

„2 n -1 



r 2 r 2 



°0 

Qi 



(27) 



(28) 



B. Fourier Transforms 



The final step of any time domain method is a Fourier transform of the time dependent information into the 
frequency domain, and for our FDTD method there are a few points which need to be made clear. Firstly, it is 
desirable to eliminate the zero frequency, longitudinal mode from our results. This is done by subtracting off the 
static part from the time dependent fields which we have calculated so that their time average is zero. Next, we must 
ensure that we obtain the properly causal solution to Maxwell's equations when we replace the integral in the Fourier 
transform with a discrete sum. We do this by adding a small positive, imaginary part 5 to the frequency. 

N t 

f(t) exp [tut] dt^Y^ f(nSt) exp + iS) (n ± l/2)5t] (29) 

n=l 

The size of this imaginary part is determined by the total time interval over which the fields are integrated. In order 
to be properly causal, the final term in the sum must tend to zero. This is of less importance when finding a band 
structure but is critical if we are to calculate the Green's function correctly. Another small detail - because of the 
choices we have made in our approximations for the time derivatives for E and H, it is necessary to include a half 
time step offset, minus for the E field and plus for the H. 



C. Conserved Quantities 



An obvious question to ask is whether our discrete Maxwell's equations have conserved quantities analogous to those 
for the continuum equations. The simplest to start with is the conservation of charge. If we consider the quantity, 

V- ■ e(r)E(r) (30) 

And then calculate the discrete version of its time derivative, 
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A- • V- • e(r)E(r) = V~ ■ A; ■ e(r)E(r) 

= v; 



V- x H' 







(31) 



Because our approximations to Maxwell's equations have preserved the form of the curls, the quantity V~ • e(r)E(r) 
remains an exactly conserved quantity even on the discrete lattice. This is the direct equivalent of V ■ D being 
conserved in the continuum case, in other words, of the conservation of charge. The same follows for V+ ■ /t(r)H'(r) 
corresponding to the conservation of the 'magnetic charge', V • B. 

The next task is to find the correct form for Poynting's theorem on the lattice. This again follows in a straightforward 
way in our formalism if we begin from the following form for the energy density. 



U(t) = - [e e(r)E*(r, t)E a (r, t - St) + HonWK ( r > 1 - 5t)H a {v, t - St)} 



(32) 



This has the correct form for the energy density in the limit St —> 0. Substituting the expressions for the generalised 
co-ordinates we obtain, 



U(t) 



U 



e afj E* a (r, t)E {v, t - St) + (-^-f ^ H'* a (r, t - St)H' p (r, t - St) 



(33) 



where Uq — Qo £ o/(QiQ2Q3l u i ■ u 2 x u 3|)- Poynting's theorem can then be obtained by considering the time differ- 
ence operator A t U(t) = (U(t + St) - U(t))/5t. This gives, 

U 



^tU(t) = ^ ( [H?(t - c,t) - H'*{v - b,t)j E 1 (r,t) 
H' 3 *(T-a,t)-H»(r-c,tj\ E 2 (r,t) 
H[*(r-b,t)-H?(r- a ,t)] E 3 (r,t) 
H' 2 (v - c,t - St) -H' 3 (r-b,t-St)] E*(r,t) 
H^(r-a,t- St) -H[(r-c,t- St)] E* 2 (r, t) 
H[(r-b,t- St) - H' 2 {t - a, t - St)] £g(r, *) 
E 2 (r + c,t)-E 3 {r + b,t)] H[*(r,t) 
E 3 (r + a,t) - E^r + c,t)\ H' 2 *(r,t) 
Ei{r + b,t) - E 2 (r + a, t)l H' 3 *(r, t) 
£ 2 *(r + c,t)-£*(r + b,i)] H[(r,t) 
E*(r + a ,t)-E*(T + c,t)\ H' 2 {v,t) 
E{ (r + b, t) - E* (r + a, t)] H' 3 (r, t) } 



(34) 



The important point to notice is that when this quantity is summed over a set of neighbouring lattice points only 
the terms associated with the surface of the integration region survive, all the volume terms cancel. This allows us 
to identify the right hand side of equation ( [34] ) as the Poynting vector integrated over the volume surrounding one 
lattice point. 



III. OBTAINING THE GREEN'S FUNCTION ON THE LATTICE 



We turn our attention now to consider how we can define the Green's function within our discrete real space 
formulism. We will begin from the continuum limit, and write Maxwell's equations as 



M 



LOP 



(35) 
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where, 



M 



+iVx 
-iVx 



e(r)e 
m("*)Mo 



We now define a six vector, 



as an eigenfunction of equation |35| with an eigenvalue u> s . We choose to normalise the F's such, 

6 6 



/u u 



And the completeness relation gives us, 



E ^(r)^" ( r ')^'i» = StfStr - r') 



We can define a Green's function in the usual way, 

s,j" 

The Green's function clearly obeys the equation, 



fK»,r)Ft ( 8 /)P. 



i<5 



(w - P- 1 M)G^(o;, r, r') = 5«*(r - r') 



We can now Fourier transform to obtain the Green's function in the time domain, 

" + 0O 



Gf r (t, r, r') = — J Gf jt (w, r, r>-*" * 



so that, 



= -iY d F s , j {T)Fl j „{T l )P j „ j ,e 

sj" 



i^g - P^M ) G R (t, r, r') = 6(t)5(r - r') 



(36) 



(37) 



(38) 



(39) 



(40) 



(41) 



Now we turn to consider how to repeat this procedure for the discrete case. We begin from equation 
the substitutions for k, uj presented in the previous section, 



x , E(r,*)\_ (A+ 
M| H(r,i) )-*[ A t - 



E(r,f) 
H(r,t) 



where M is now 



M 



+iV~x 
-iV _ x 



(42) 
(43) 

3q, but apply 

(44) 

(45) 



We want to construct the G R (t) which will obey the equation 



G fl (t,r,r') = W(r-r') 



(46) 



Try constructing, 
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G? jt (t,T,T>) = -i«yt^F„-(r)^ // (r')P J -„ i , e -^(*-«);i > 
G^,(t,r,r')=0 ;i<0 



then, 



G* (u,,r,r') = £g* (t.r.iV* 

t=<5t 



(47) 
(48) 



1 _ e i{u-u s )N t 5t 



If N t is large and cj = u + iS, 



Gf } , (u, r, r') - -tft ]T F sJ (r)^.„ (r'jP,-, 
s,i" 



■iuj s 8t 



-i(uj— u) s )8t ^ 



(49) 



which in the limit St — > gives the correct form for the frequency dependent Green's function. So simply by setting 
t = St in equation [l?] to give, 

Gfj, (t = St, r, r') = -i St S jf S(r - r') (50) 

as the appropriate starting condition and applying equation (^) we can calculate G R {t) at all subsequent times. 

Once we have found the Green's function it is a simple matter to calculate useful physical quantities from it, such as 
the photonic density of states. The local density of states, for example, is found in the usual way, from the imaginary 
part of the trace of G fl (w)Ej. 



-Im 



w.r.r 



(51) 



Similarly, the band structure can be easily determined by locating the poles in the Green's function and so identifying 
the normal modes of the system. 



IV. RESULTS 



We shall now put the ideas of the previous sections to work for two different physical systems. 



A. Bragg Stack 

In the first case we will look at what is probably the simplest one dimensional photonic crystal that can be imagined, 
the dielectric multilayer, or Bragg stack. This crystal is formed by stacking together alternating layers of dielectric 
of high and low refractive index (see figure |l|). Each layer is of infinite extent in the plane and we choose parameters 
such that the high refractive index (n = 3.6) planes have a thickness of 0.3a and the low index (n = 1.0) planes 
have a thickness of 0.7a where a is the lattice spacing. This choice leads to a sizable band gaps for electromagnetic 
waves propagating normal to the planes as can seen in figure ^[ The frequencies in this figure are scaled so as 
to be dimensionless, what we actually write for the frequency is wa/(27rco). The wavevectors are also written in 
dimensionless units of k.a so that the edge of the first Brillouin zone occurs at ±tt. 

In figure ^| we show the local density of states for a point inside the Bragg stack. Again the frequencies are given in 
dimensionless units and the density of states itself is in arbitrary units. The band gaps are clearly visible as frequency 
windows over which the density of states is strongly suppressed. At the band edges the Van Hove singularities are 
also apparent corresponding to the points on the band structure at the band edges where du/dk tends to zero. 

We can go one step further and add a defect layer to the otherwise perfect crystal. We do this by introducing a super- 
cell consisting of 25 repeat units of the Bragg stack, then a defect of 0.3a of the high refractive index material, and 
then another 25 unit cells. The local density of states at the centre of the defect can be seen in figure [| superimposed 
over the density of states for the perfect crystal. The dominant feature are the new peaks which have appeared in the 
band gaps. These correspond to localised modes associated with the defect. That these modes are tightly localised 
can be easily shown by calculating the local density of states in the crystal several lattice spacings away from the 
defect. The perfect crystal result is recovered. 
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B. Diamond Lattice 



The second system which we shall consider is a fully three dimensional photonic crystal made up from air spheres 
embedded in a dielectric host arranged in a^diamond lattice structure. This structure was one of the first to be 
shown to have a complete photonic band gapQ. Figure || gives a diagram of the diamond structure and indicates the 
unit cell which we choose. This unit cell is not primitive but is chosen for convenience with lattice vectors; a = ai, 
b = (a/2)i + (\/3a/2)j and c = \/6ak where a is the distance between nearest neighbour lattice points. The exact 
parameters we have chosen lead to a fairly sizable band gap. The refractive index for the host material is n — 3.6, 
the spheres have n = 1.0, and the ratio of sphere radius to the distance separating nearest neighbour spheres is 0.43. 
The band structure for this system in the directions T — K and T — X is given in figure |^ where it is clear that a band 
gap is opening up over a frequency range of about 3.2 and 3.9 in dimensionless units. This band gap is confirmed in 
the density of states, in figure [?]. The density of states calculation is performed by creating a super-cell by stacking 
unit cells together, eight in each direction. Because each of our unit cells is not in fact primitive but itself contains 
three primitive cells, the super-cell consists of 1536 primitive units cells. By creating a larger unit cell we reduce the 
number of k-points which we need to sample in order to build up the complete local density of states. 

Next we can introduce a defect into the diamond lattice by again creating a super-cell by stacking unit cells 
together and then cutting a hole in one of the dielectric parts of the cell. This defect will have a localised defect mode 
associated with it and by altering the amount of dielectric material we cut away we can tune the frequency of the 
mode. We choose the size of the defect so that the defect mode lies within the band gap for the crystal. We cut away 
a parallelepiped region spanning adjacent unit cells starting from (a/2, b/3, c/3) and with edges of 5a/6, 5b/6 and 
5c/12. Figure || shows the density of states for the super-cell with the defect and as for the Bragg stack case, several 
defect modes are clearly seen. We have also checked that we we have made the super-cell sufficiently large by testing 
that the position of the defect mode peaks do not disperse with k. 



V. CONCLUSIONS 



In this paper we have shown how an Order-N, finite difference time domain method can be used to calculate the 
Green's function in a simple and straightforward way. Specifically, our new formalism enables us to give a simple 
derivation for the numerical stability criterion, exact statements of charge and energy conservation and allows us to 
use non-orthogonal co-ordinate systems without the usual computational overheads. 

From the Green's function a whole range of other physical quantities can be found such as the local density of 
states. The times required to calculate these quantities scales linearly with the size of the system, so our technique 
is of particular importance in analysing systems with very large or complicated unit cells. Specifically, we have 
calculated the local density of states for both a one and three dimensional photonic crystal containing a defect and 
have recovered the expected result. The defect has the effect of introducing a highly localised mode the frequency of 
which is determined by the size of the defect. By a careful choice we have found defect states within the photonic 
band gap for the crystal. These localised states have potential application in photonic cavity laser structures where 
the efficiency of the laser is enhanced by suppression of emission into all but the lasing mode. The linear scaling 
of this method with system size should allow modelling of realistic designs for cavity structures without prohibitive 
computational overheads associated with traditional computational schemes. 

The computer program used to calculate the results presented in this paper has recently been submitted to the 
CPC International Program Libraryll3. 



1 K. S. Yee, IEEE Trans. Antennas Propag. 14, 302 (1966). 

2 C. T. Chan, Q. L. Yu, and K. M. Ho, Phys. Rev. B 51, 16635 (1995). 

3 J. B. Pendry, J. Phys.: Condens. Matter 8, 1085 (1996). 

4 E. Yablonovitch, J. Phys.: Condens. Matter 5, 2443 (1993). 

5 E. Yablonovitch and T. J. Gmitter, Phys. Rev. Lett. 63, 1950 (1989). 

6 H. Hirayama, T. Hamano, and Y. Aoyagi, RIKEN Super Computing Prog. Rep. 1, 1 (1996). 

7 K. M. Leung and Y. F. Liu, Phys. Rev. Lett. 65, 2646 (1990). 

8 Z. Zhang and S. Satpathy, Phys. Rev. Lett. 65, 2650 (1990). 

9 K. M. Ho, C. T. Chan, and C. M. Soukoulis, Phys. Rev. Lett. 65, 3152 (1990). 



9 



10 J. B. Pendry and A. MacKinnon, Phys. Rev. Lett. 69, 2772 (1992). 

11 J. B. Pendry, J. Mod. Optics 41, 209 (1994). 

12 F. Wijnands et al, Optical and Quantum Electronics 29, 199 (1997). 

13 J. Lee, R. Palandech, and R. Mittra, IEEE Trans. Microwave Theory Tech. 40, 346 (1992). 

14 A. J. Ward and J. B. Pendry, J. Mod. Optics 43, 773 (1995). 

lj E. N. Economou, Green's Functions in Quantum Physics ( Springer- Verlag, Berlin, 1990). 

16 A. J. Ward and J. B. Pendry, A Program for Calculating Photonic Band Structures and Green's Functions using a non- 
orthogonal FDTD Method, Submitted to Computer Physics Commun. 

FIG. 1. A dielectric multilayer or Bragg stack formed by alternating layers of high (shaded) and low (unshaded) permittivity 
materials. The thickness of the high dielectric layer is d and the lattice period is a 



FIG. 2. The photonic band structure for Bragg stack 



FIG. 3. The density of states for the ideal Bragg stack. The Van Hove singularities in the perfect crystal are clearly seen 

FIG. 4. The density of states both with and without a defect inserted into the ideal Bragg stack. Both the Van Hove 
singularities in the perfect crystal and the localised state associated with the defect are clearly seen 

FIG. 5. Diagram of the diamond lattice. The bravais lattice is FCC and the lattice points are shown as the black dots. On 
each lattice point are placed two 'atoms' - dielectric spheres in our case; one at (0, 0, 0), one at (0, 0, c/4), where c = \/6a. Some 
of the 'atoms' are shown as the grey circles. The unit cell we model is the parallelepiped enclosed by the thick black lines 

FIG. 6. Partial band structure for spheres on the diamond lattice in the F — K and F — X directions. The presence of the 
band gap is clearly shown 



FIG. 7. The density of states for the diamond structure 

FIG. 8. The density of states both with and without a defect for the diamond structure. The localised state associated with 
the defect is clearly seen 
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